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The decay of unstable states when several metastable states are available for occupation is 
^3 \ investigated using path-integral techniques. Specifically, a method is described which allows the 

^^ . probabilities with which the metastable states are occupied to be calculated by finding optimal 

paths, and fluctuations about them, in the weak noise limit. The method is illustrated on a system 
described by two coupled Langevin equations, which are found in the study of instabilities in fluid 
dynamics and superconductivity. The problem involves a subtle interplay between non-linearities 
and noise, and a naive approximation scheme which does not take this into account is shown to 
be unsatisfactory. The use of optimal paths is briefly reviewed and then applied to finding the 
conditional probability of ending up in one of the metastable states, having begun in the unstable 
state. There are several aspects of the calculation which distinguish it from most others involving 
r^ ■ optimal paths: (i) the paths do not begin and end on an attractor, and moreover, the final point 

O ' is to a large extent arbitrary, (ii) the interplay between the fluctuations and the leading order 

•^ , contribution are at the heart of the method, and (iii) the final result involves quantities which are 

not exponentially small in the noise strength. This final result, which gives the probability of a 
particular state being selected in terms of the parameters of the dynamics, is remarkably simple 
C^ , and agrees well with the results of numerical simulations. The method should be applicable to 

ri ' similar problems in a number of other areas such as state selection in lasers, activationless chemical 

reactions and population dynamics in fluctuating environments. 
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I. INTRODUCTION 



K*" ' The decay of metastable states, due to thermal or other random fluctuations, is a phenomenon seen in many diverse 



areas of science, and consequently has a huge literature associated with it. In the simplest cases, where a potential 
J;^ ■ can be defined and the states assigned a particular value of this potential, the decay process can be viewed as noise 
I — ' activation over the potential barriers that separate the metastable state under consideration from all of the other 
f^ ' accessible metastable states of the system. The average time taken to escape from a potential well (i.e. for the state 
T— I ' to decay) is of the order of exp /SV/D, where /S.V is the height of the barrier to be surmounted and D is the strength 
^^ : of the noise. Thus the picture we have in this case, is of a set of metastable states, with transitions between them 
.4.^ which occur with probabilities which depend on the nature of the potential between the states and on the strength of 
3 . the noise. 

By contrast, the decay of unstable states, although a similarly widespread phenomenon, has been studied much 
less. In terms of the above picture, the system starts at or near a maximum of the potential and makes transitions 
to the accessible metastable states of the system with various probabilities, which as before depend on the nature of 
Q \ the potential between the unstable and metastable states and on the noise. In an earlier paper M we introduced a 
Q . scheme which enabled us to calculate the probabilities with which the various metastable states are selected. Our 
aim here is to extend this work, by giving a fuller presentation of the ideas and techniques involved, justifying some 
of the earlier approximations which were made and discussing the link with real systems in more detail. 

The phenomenon of the selection of metastable states, from states which have become unstable, is at the heart of 
pattern formation and the origin of complex structures, with the selection of non-trivial states being governed partly 
by the deterministic dynamics and partly by the noise acting on the system. Thus the picture painted above, while a 
simplified version of this general scenario, contains many of its essential features. In fact, the above structure can be 
derived from the equations describing the entire system by focusing on the unstable modes and the modes which have 
the potential to be selected, and treating all of the other modes as background noise. The resulting dynamics then 
comprises a small number of coupled ordinary differential equations acted on by noise. If this flow is potential, then 
the above picture is recovered; if not, our techniques are still applicable, but it becomes more difficult to visualize. 

State selection of the kind we have been describing is ubiquitous. In fluid dynamics, it appears when Rayleigh- 
Benard convection rolls of given wavenumber are formed after the decay of unstable ones. A set of equations describing 



the non-linear coupling between parallel rolls of different wave numbers may be derived B , and when the effects of 
the other modes are incorporated |3| , a set of coupled stochastic differential equations of the kind mentioned above 
are generated. In superconductivity, exactly the same equations as in the above example govern state selection 
in, for example, narrow superconducting rings. This is simply because the amplitude equation which governs the 
instabilities in the fluid dynamics example is nothing but the Ginzburg-Landau equation for a superconductor Wj . A 
mode truncation then gives exactly the same set of coupled differential equations acted upon by noise . In chemical 
kinetics, it has become clear in the last decade or two, that there are many important chemical reactions in which 
a barrier to the formation of an excited state is not present. Examples of these activationless reactions |g] are the 
electronic relaxation of triphenylmethane dyes and barrierless electron transfer in solution. The potential discussed 
above is a reaction potential energy surface in this case 0. In lasers, such as the ring dye laser [||, decay of an 
unstable mode to metastable ones can occur under the right operating conditions. The coordinates in this case are 
the mode amplitudes and Langevin-type equations are derived within the semiclassical theory of the laser H. In 
population dynamics, the Gause model of two competing species [H again falls into this class. When the competition 
is played out with a fluctuating environment, the resulting stochastic differential equations once again fall into the 
generic class that we have been discussing |1^. It should, however, be noted that in this case there is no potential 
and moreover the noise is multiplicative. 

Although many of these phenomena have been known for some time, the investigation of the state selection aspect 
has been hampered by the lack of a suitable calculational tool. In the case where a potential V exists, the noise is 
additive and only two modes, x and y, are considered, the equations take the form 



dV 

dx 
dV 



Vxit) 



y = —g^+Vy{t), (1) 

where r]x and rjy are white noises of strength D. But, even in this, the simplest non-trivial case, these equations are 
difficult to study mainly because, as we shall discuss later, in the region where state selection occurs, the coupled nature 
of the equations, their non-linearity and the noise, are all important. However, as we have shown ||I|, there is a method 
that can take all these aspects into account, and that is the path-integral formulation of stochastic dynamics. This 
method succeeds where others fail because the equations can be represented as a path-integral without approximation] 
systematic approximation techniques for path-integrals developed over the years can then be used as the basis of a 
calculational scheme. 

Most of the previous theoretical work on this problem has been limited to systems with one degree of freedom, 
governed by the single equation x — —V'{x) +r]{t), because many of the complexities mentioned above are not present 
in this one-dimensional case. Suzuki and co-workers developed a theory for the decay of an unstable state in one 
dimension in a series of papers fl^, as did a number of other authors [p"3|-p^. However, the one-dimensional theory, 
although much easier to deal with, has none of the subtleties inherent in state-selection in higher dimensions: the 
decay is simply either to the right or to the left. Some studies |20|-p2[ purported to go beyond one dimension, but 
in fact considered spherically symmetric potentials, so that the problem could be reduced to a quasi-one-dimensional 
problem in the radial coordinate. Once again, the resultant structure is not rich enough to address the question of state 
selection. Other approaches, such as replacing the stochastic differential equations by the corresponding deterministic 



equations, but with random initial conditions |23|, are unsatisfactory in other ways. Probably the investigation nearest 
to our own has been by Mangel [e3, however his main interest was not state selection. 

The plan of the paper is as follows. In Sec. II we discuss the assumptions underlying the picture of state selection 
which we have outlined and introduce a generic model which we use to describe our calculational scheme in detail. 
An important aspect of our method is the realization that, typically, a state is selected well before the system reaches 
this chosen state. This fact is used to simplify the problem in Sec. HI, where it is also shown that naive calculational 
prescriptions, based on linearization of the initial dynamics, fail. A more systematic approach based on path-integral 
techniques is introduced in Sec. IV and the relevant optimal path is determined. This gives the leading order 
contribution in the limit D —f 0. The next to leading order contributions are determined in Sec. V. The results of this 
analytic approach are given and compared to Monte Carlo simulations in Sec. VI and our conclusions are presented 
in Sec. VII. There are two appendixes. Appendix A describes some of the elementary approaches described in Sec. 
HI in more detail. Appendix B contains technical aspects relating to the determination of the y-coordinate of the 
optimal path and the calculation of the action of that path. 



II. GENERAL CONCEPTS 

We have already given several examples of situations where state selection, following the decay of an unstable state, 
occurs. In this section we will explore the different ways in which the initial state of the system could arise, that is, 
the origin of the unstable state, and give an intuitive description of the decay and subsequent state selection, which 
will form the basis of our analytical approach. 

An unstable state may arise by mechanisms which are either non-adiabatic or adiabatic. In the non-adiabatic case, 
the system, in an initially (meta)stable state, is transported very rapidly to an unstable state. The time scale for 
the transition from the stable to unstable state is much more rapid than the characteristic time scale defining the 
natural dynamics of the system. In this context the transition from the stable to the unstable state can be ignored 
and we simply characterize the system as having been prepared in an unstable state. Examples of this type were 
mentioned in Sec. I and include chemical reactions with no activation barrier and population dynamics in a fluctuating 
environment. Quasi-one dimensional superconductors |5[| are another example. In the case of the chemical reactions, 
an ultra-fast light pulse (a femtosecond laser) pumps molecules into an excited state. In terms of potential surfaces, 
the molecules are initialized in an unstable state on the reactive surface. The subsequent reaction can be viewed as 
a nuclear rearrangement on this reactive potential energy surface. The evolution of the reaction, that is the relative 
preponderance of reactants and products, can be studied using other ultrafast lasers. In the population dynamics 
example, the population is initialized so as to consist of only a few individuals of each species. During the initial 
period the population of both species grow exponentially, independently of the other, until competition drives the 
population of one of the species to zero. 

In the adiabatic case, an unstable state may arise when a (meta)stable state is transported through the instability 
over a time scale that is slow compared to the time scale characterizing the natural dynamics of the system. An example 



of this process is provided by a quasi one-dimensional superconductor as studied in Ref. |25|. In this case, the system is 
driven by a voltage source that accelerates the supercurrent. However, this acceleration cannot continue indefinitely 
and the system becomes unstable, a situation associated with the critical current of the superconductor. Once 
unstable, the relaxation process occurs via "phase slips" in which spatially localized regions of the wire temporarily 
lose their superconducting properties and carry "normal", i.e. non-superconducting, current. This process dissipates 
the excess energy and the system relaxes back to a locally stable state. Mesoscopic wires, i.e. wires that are not in the 
thermodynamic limit, have a finite number of discrete metastable stationary current carrying states. State selection in 
this case is characterized by the competition between these metastable current carrying states. The superconductivity 
instability described above is an example of the Eckhaus instability. 

Eckhaus instabilities arise in many physical systems in addition to quasi one-dimensional superconductors, including 
fluids, nematic liquid crystals and lasers. In these systems, stationary one-dimensional periodic patterns are stable 
for a range of wave vectors, Q. For instance, in the case of the common Eckhaus instability, stationary solutions exist 
for Q^ < 1, but are only stable if Q^ < 1/3. By changing the control parameter slightly, Q-states with Q^ < 1/3 may 
be shifted into the unstable Q^ > 1/3 regime. The essential features of the subsequent changes which follow from 
this may be understood by performing a mode-decomposition of the relevant amplitude equation and keeping only 
the previously stable mode and the destabilizing modes with the largest growth rates S.An adiabatic elimination of 
the unstable mode then leaves us with coupled ordinary differential equations for the amplitudes of the destabilizing 
modes. This adiabatic assumption is equivalent to assuming that the form of these differential equations does not 
change on the time-scale of state selection, i.e. in the time taken for one of the destabilizing modes to outcompete 
the others |g^ 

An obvious question which now arises is the following: how do we know that the system is initialized in a state where 
the amplitudes of the competing modes, final products, etc are zero? Or expressed in the topographical terminology 
of Sec. I: How do we know that we begin exactly at the top of the hill? Here we are assuming that the small non-zero 
amplitude, required to initiate the growth, is induced by fluctuations i.e. by the noise. Once again, there are several 
possibilities. It may be that in every realization of the stochastic process the system starts near, but not necessarily 
at the origin. As long as there is no bias favoring any one of the competing modes, an average over many realizations 
of the process will give the same results as if the system had started at the origin in each case. Alternatively, the 
initial probability distribution may be so sharply peaked about the origin, that any small initial bias is irrelevant to 
the ultimate choice of state made by the system. Yet another possibility is that it is just not possible to start exactly 
at the origin (the example taken from population dynamics is a illustration of this). In these cases the sensitivity of 
the final result to changes in the initial starting point will have to be carefully investigated. In addition to all of these 
specific reasons, on general grounds beginning at the origin seems very natural, since it simply means that none of 
the final products or final modes are initially present. 

Since there appear to be good reasons to initialize the system at the origin in many cases, it will be assumed in the 
subsequent theoretical development. In all of the examples given in Sec. I, the competing modes grow exponentially 



for an initial period, without any significant interaction with each other. This corresponds to a hnear growth law of 
the form xi = aixi, where xi is the amplitude and ai the growth rate of the Ith mode. Clearly, the main interest 
from a state-selection point of view, is in the nature of the non-linear interaction terms between the modes. Since 
the purpose of this paper is to give a clear and intuitive introduction to our calculational scheme, we shall choose to 
illustrate our method on an example where only two competing modes are present, and where the interaction terms 
may be derived from a potential. In other words, the governing equations will have the form (|l|). Obviously, two is 
the minimum number of modes we need in order to discuss state selection, but having only two modes moving on a 
potential surface has the advantage that we can describe our method using the pictorial language of "hills and valleys" 
which we have already adopted on several occasions. Furthermore, since as we have already mentioned, the equations 
derived in both the fluid dynamics and superconductivity examples are identical |g,^ we will work with these. They 
are 

X — ax — "fxy^ — 5x^ + rjx{t) (2) 

y = /3y - ^yx'^ - ey^ + riy{t) , (3) 

where x and y are the two competing modes with growth rates a and /3 respectively. All five parameters characterizing 
the model (a,/3, 7,(5 and e) are assumed to be positive. The noise terms rix{t) and tjyit) are taken to be Gaussian 
random variables with zero mean and 

{rj,{t)Tj,{t'))^2DS{t-t'), (4) 

where i and j are either x or y and D is the noise strength. Equations (0) and (13) fall into the class of those given by 
(|l[), since they are derivable from the potential 

V{x,y)^-^x'-^y' + lxV + l^' + ly'. (5) 

We wish, however, to stress that the method is not restricted to systems with only two modes nor is the existence of 
a potential a prerequisite. This will become clear as the method is explored in subsequent sections. 

It is interesting to note that the simple generalization of the model in which the strengths of the noises rji (t) and 
r]2{t) are both of the same order, but not equal, already leads to a more complicated situation. If we suppose that the 
strengths of these noises are Di and D2 respectively, then we may introduce Ci{t),i = 1,2 such that the correlation 

function of these new noise terms is exactly (0), by writing rii{t) = p^ Ci{t) where the pi are constants given by 
Di = piD. This rescaling of the noise may, in turn, be eliminated from Eqns. (^ and (^ by rescaling x(t) and 
y{t) by defining new variables xi and X2 via x{t) — p^ x\{t) and j/(i) — p.{ X2{t). New constants replacing S and e 
which absorb this change of scale can easily be defined, but the interaction terms jxy^ and 'yyx'^ become 7i2a;ia;2 and 
721X20;^ respectively, with 712 ^ 721- Equations such as this are not derivable from a potential, nevertheless minor 
modications of our method are applicable to this case. 

The potential V{x, y) for a particular choice of parameters is shown in Fig 0(a). 
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FIG. 1. (a) Illustrative plot of V{x,y) with a — j3 — 2,^ = 4, 5 — e= 1/5. (b) Contours of zero force. 



Before discussing the dynamics of the system, we have to identify the extrema of the potential. These are given 
by the intersection of the two sets of curves dV/dx — and dV/dy = or, exphcitly, by the intersection of 
{x = 0, Sx^ + -fy"^ = a} and ly = 0, -fx'^ + ey^ ~ 13} . These two elhpses, along with the x and y-axes, are shown in 
Fig |](b). By comparing Figs |l|(a) and Fig |l|(b) we see that there is (i) a maximum at the origin (denoted by a filled 
circle), (ii) four minima at the intersection of the ellipses and the axes (denoted by open circles), and (iii) four saddle 
points at the intersection of the two ellipses (denoted by crosses). The precise positions of these extrema are: 



Minima: (x, y) = ("i (^)'^' , o") ; ix,y)^U ±(^) j. (6) 

Saddles : (x, y) - J-^ (±^Pj-ae, ±^aj-ps) . (7) 

^7-^ — de ^ ^ 

In fact, there are some conditions which need to be imposed on the parameters of the model if saddle points are to 
exist. These are that 137 — (US, (Uj — ae and 7^ — Se should all have the same sign. Since we will normally assume that 
the stability parameters, S and e, which ensure that the potential is bounded below, are small, we will take Pj > ae 
and 017 > pS, which implies that 7^ > Se. 

The dynamics given by (|l|), (j^ and (|^) is equivalent to an overdamped particle moving on the potential surface 
shown in Fig [^(a) and which is also acted upon by white noise. In a typical realization, the particle will begin at 
the maximum and perform a random walk which is in the vicinity of the origin at early times, but which eventually 
explores an ever larger region. Eventually, the particle gets far enough from the origin that the deterministic dynamics, 
specified by V, begins to have a significant impact and the particle accelerates down towards the saddles and the 
minima. In this paper our main concern is state selection: we are not concerned with the approach the particle makes 
to a particular minimum, since by this stage this minimum will almost certainly be the selected state. For this not 
to be so, the particle would have to hop over a barrier — an extraordinarily rare event. In fact, as is clear from Fig 
0(a), as soon as the particle has passed the x or y coordinate of one of the saddle points, it has effectively chosen the 
final state and its subsequent motion is of little interest to us. We therefore arrive at the conclusion that the saddle 
points are the major factor influencing state selection and by comparison the minima are of little consequence. This 
can be made more transparent if we imagine that S and e are very small, so that the minima m) are now at very large 
values of x and y. By contrast, the positions of the saddle s ha ve cha nged much less: in the limit where S and e go to 
zero their x and y coordinates tend to the finite values \fWhl ^^^ 1/0^/7, whereas the minima tend to infinity. The 
topography now consists of long valleys, leading from the head of the valley between two saddles, all the way down to 
a minima at the end of the valley. As soon as the particle enters the valley it is extremely unlikely to escape over the 
sides, and so extremely unlikely it will do anything else but move to the minima at the end of that particular valley. 
This makes it very clear that it is the saddles, and not the minima, that we should focus on if we wish to understand 
state selection. 

These ideas are well illustrated by plotting out a few Monte Carlo trajectories. Such simulations are extremely 
simple to carry out, and later we will compare our analytic expression for the probability of a particular state being 
selected with the proportion of runs which ended up in that state. For the moment, however, we are interested in 
the nature of individual trajectories. The five which are shown in Fig 2 are typical: the particle carries out Brownian 
motion about the origin, to a greater or lesser extent, and then selects a particular valley. Even after state selection, 
there may still be large deviations, but eventually the particle settles down near to the axis along which the valley 
runs. 




X 



FIG. 2. Typical trajectories for the full potential V{x,y) with a = l3 — "/ = 1, 5 = e = 0.1, and with noise strength D — 0.01. 

III. THE REDUCED POTENTIAL 

We have seen in the last section that the essential features of state selection become much clearer in the limit 
(5, e — *■ 0, when the saddle points are completely separated from the minima, which have moved off to infinity. For the 
rest of the paper we work in this limit in order to illustrate our method in the simplest possible way. We will call the 
problem defined in Sec. II when the limit is taken, the reduced problem. The reduced potential is then defined to be 



VR{x,y) = -2^^ 



/3 2 , 7 2 2 

2^ +ry 



(8) 



A plot of Vr looks very similar to the plot of the full potential V in Fig H(a), since this figure emphasizes the little- 
changed central features of the potential. The main difference is that whereas in V the valley floors start to ascend 
again at large |a;| and |y|, in Vr they keep descending. However, as was made clear in Sec II, this behavior is of no 
interest to us, and the fact that V and Vr differ in this way is immaterial. 

From (0) , we see that Vr has only two types of extrema: (i) a maximum at the origin, and (ii) four saddle points 
at {x,y) = (±X,nin,±Fmin) whcre 




7 



(9) 



The reason for the subscript "min" is that the values given by (g) are the smallest values of |a::| and \y\ at which we 
can assume that state selection has taken place. A simple stability analysis near these saddles shows that the stable 
and unstable directions are at an angle of 7r/4 to the axes. For example, in the positive quadrant, if the particle 
approaches the saddle along the line of steepest descent y = x, it will have a tendency to move down the slopes which 
arc to the right and left, rather than carry on through the saddle and straight up the incline ahead. 

Just as the plot of Vr is similar to that of V if we cut it off near to the head of the valleys, so the Monte Carlo 
simulations are similar in both cases if we do not follow the trajectories too deeply into the valleys. In Fig. 3, the range 
of X and y used to display the trajectories are smaller than those in Fig. 2, so that the central region is emphasized. 
This makes the initial Brownian dynamics seem more obvious, but in reality the nature of the trajectories before and 
during state selection are essentially identical to those of the full problem shown in Fig. 2. 
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FIG. 3. Typical trajectories for the reduced potential Vr{x, y) with q = /3 = 7 = 1, and with noise strength D = 0.01. 



In the next section, we will develop a calculational scheme, based on the path-integral formulation of the stochastic 
process, to determine the probabilities of entering the x- or y-valleys as a function of the model parameters a, /3, 7 
and D. However, one might feel that such a sophisticated theory is not necessary: it should be possible to obtain 
a satisfactory theory by constructing an approximation based on the picture we have built up in this and the pre- 
vious section. We will therefore end this section by constructing an example of such a theory and show that it is 
unsatisfactory for a variety of reasons. 

A simple theory of state selection might have the following ingredients: 

1. In the initial period the noise and linear growth of the modes dominate. Therefore, neglect the non-linear 



interaction between the modes (i.e. set 7 = 0). 



2. State selection will be specified in the following way. Define four sectors in the xy- plane by drawing lines through 
the origin and the saddle points. The particle will select the state corresponding to a particular sector if it is in 
that sector for x « Xram and y w ymin. 

A calculation based on these two assumptions is given in Appendix A. Even within the framework that these 
provide, there turn out to be many possible variants, each giving slightly different results. Leaving this aside for the 
moment, it is shown in Appen dix A that in one of the simplest schemes along these lines, the probability of ending 
up in an x- valley is (see Eqn. ( A14 )) 



N, 



1 + DP 



(10) 



where 



~ -/D [a~(3) 
D = — - and p = 

ajj a 



(11) 



This has the correct qualitative behavior: when a — 13. 
unity. A more careful calculation gives (Eqn. (A15)) 



that is, p = 0, Nx = 1/2. As p increases, N.j. increases towards 



N^^ - tan^i b-P . 

IT 



(12) 



which again has the correct qualitative features. 

We now compare these two formulas with the results of Monte Carlo simulations. We took P,j and D to be fixed 
in the simulations, at values 1.0, 1.0 and 0.01 respectively, and varied a from 1.0 to 3.0 in steps of 0.2. For each value 
of a we performed a large number of runs where the particle started at the origin and ended up in one or other of the 
valleys (for more details see Sec. VI). We then simply counted the number of times that the particle ended up in a 
particular valley and expressed this as a fraction of the total number of runs made. The results are shown in Fig. 4. 
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FIG. 4. Probability of flowing into an a;-valley as a function of a, with /? = 1, 7 = 1 and D = 0.01. The dots show simulation 
results, the full curve is the result (FuJ) and the dashed curve is result (p^. 

It is clear that the results of the naive approach based on assumptions 1 and 2 above, worked out in Appendix A 
and given by (nO) and (|l3) are not quantitatively correct. They mainly overestimate the probability of the particle 
ending up in the a;- valley. In addition, Eqn. (p^, which represents a more refined calculation, actually compares less 
favorably with the simulation results, than does (^0|). This clearly points to the unsatisfactory nature of this simple- 
minded scheme, as one would expect the more refined calculation to compare more favorably with the simulation 
results. 

We do not wish to pursue techniques based on this approach further in this paper. It was included specifically to show 
that there is a delicate interplay between the non-linearity and the noise in (g) and (|3|) and any attempt to incorporate 
one or the other of these using an unsystematic approximation scheme, as above, is likely to give disappointing results. 
We do not rule out being able to construct a particular scheme of this type which will give reasonable results: there 
are enough possible variants that it may be possible to interpolate between these by introducing free parameters which 
can then be fitted to the Monte Carlo data. However, such an ad hoc scheme is unnecessary, since we will show in the 
remainder of the paper that a systematic approach exists which gives simple formulas which are in good agreement 
with the Monte Carlo data. The method is based on the use of optimal paths, and so we begin with a brief review of 
the path- integral formulation of Langevin dynamics. 
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IV. OPTIMAL PATHS 

In this section we review the formulation of stochastic differential equations, such as (P, as functional integrals and 
obtain the dominant contribution to the conditional probability we wish to determine in the limit where the noise 
strength tends to zero. 

The conditional probability that the system is in the state (a;/,y/) at time T, given it was initially in the state 
(0, 0) at t = is 



P[xf,yf,T\0, 0, 0) = {5 (xf - x{T)) 5 {vf - y{T))) 



ic 



(13) 



where IC denotes the initial condition a;(0) = 0, y(0) = on the stochastic process and x{T) and y{T) are solutions 
of (0) and (0). The average in (O) is over Gaussian white noises rix{t) and riy(t) with zero mean and correlation 
function given by (P. In terms of functional integrals (|lj) equals 



C J^^ D^,D7^y 5 (xf - x^{T)) 5 ivf - y^iT)) exp | -^ ^ dt [r^l{t) + r;2(t)] 



(14) 



where C is a normalization constant and the subscript rj on x{T) and y(T) is to emphasize that they depend on r](T) 
through (0) and (0). Using these equations to perform a functional change of variable from {rix,riy) to {x,y) yields 



cJ^^DxDyJSixf~xiT))Siyf-yiT))e^p\-^J^ dt 
where J is the Jacobian of the transformation. Expressed as a path-integral 

r(T)=rf 



dV 
dx 



dV 
dy 



P(f/,T|d,0) =C I DrJ[r\ eiii^p {- S[r\ / D} , 

Jf{0)=0 

where r = {x,y). The action S[r\ and the Jacobian J[r] are functionals which are given by 



S[x,y] = ^ / dt 



dV 
dx 



dV 
dy 



(15) 



(16) 



(17) 



and 



J[x,y] = Det 



6ff 
Sf 



1 

oc exp < — 



dt 



d^V d^V 



dx"^ dy^ 



(18) 



For D ^ 0, the path-integral dll) is dominated by solutions of the Euler-Lagrange equations SS/6r{t) — 0, which 
satisfy the boundary conditions r(0) = and f{T) — rf. Let the solution of least action be denoted by rdt)- Then 
writing r{t) — fc{t) + Sr{t), we have 

P{rf,T\6,0) =C exp{~S{fc)/D} [ DSr J[fc + Sr\ 



xexpi-— /" dt' f dt"5r{t') 



1 



5^S 



2 5r{t')5r{t") 
Scaling 5r by D^''^ and performing the Gaussian functional integral yields 

P(r>,T|0,0) =C' exp{-5(r,)/i?} J(r,) Det ^ ^ 



5r{t") + 0{5rf 



5r{t')5r{t") 



-1/2 



[1 + 0{D)] 



(19) 



(20) 



The new overall constant, C, is immaterial since, as we have seen in Sec. Ill, we are interested only in the probability 
of ending up in the a;- valley or y- valley, and we normalize these probabilities according to Nx + Ny = 1. We will 
therefore omit the overall constant from now on. 
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While all the derivations in this section have been carried out for potential systems with two degrees of freedom, 
it should be clear that they generalize in an obvious way to systems of more than two degrees of freedom and those 
where no potential exists |2^ . We may summarize the result of performing the functional steepest descent on (|l^) to 
next-to-leading order by 



P 



(f;,r|0,0) = P«(r>,r)exp{-p(")(r>,r)/I?} [1 + 0{D)] 



(21) 



where the leading order contribution, P^^\ is just the action of the optimal path rc{t) and the next-to- leading order 
contribution is 



P(i)(f/,r) = J{r,) (Deti(r-))-i/2 



(22) 



Here L is the matrix formed from the second order functional derivative of the action functional evaluated at the 
optimal path: 



Wc) 



s^s 



5r{t')5f{t" 



(23) 



The result ( pT| ) is the starting point for our method: if we can determine the functions P^"^ and P^^\ then we will 
have a form for the conditional probability valid when the noise is weak. This in turn will enable us to obtain a 
formula for the probability that cither the x- valley or y- valley is selected, as a function of a, /3, 7 and D. We shall 
devote the rest of this section to the determination of P^^"^ and leave the calculation of P'^^^ to the next section. 

The case of interest to us in this paper is when V is the reduced potential (||) . The Eulcr-Lagrange equations for 
this problem are 



X — X [a — 72/^) — 2^xy^ (/? — 72:^) and 
y = y{P- ix^) - 2jyx^ (a - jy"^) . 



(24) 
(25) 



It is important to realize that there are two distinct dynamics associated with the problem under consideration. The 
first is the stochastic dynamics given by (|l|) with the potential (||). This was our starting point, the basis of the 
intuitive discussion of the dynamics given in Sees. II and III, and the dynamics of the Monte Carlo simulation. 
The second dynamics is the deterministic dynamics, given by ( p4| ) and (Eq), which describes the Z? — > limit of the 
stochastic dynamics. They are quite different and it is important not to carry over intuition from one to the other 
without careful consideration. From (|1^ 



SM 



dt 



\--2 

—r 

2 



U(r) 



dt 



dV 



where 



"^^^-m-m:' 



(26) 



(27) 



Since the last term in ( |2q ) is a constant, and consequently gives zero variation, the Eulcr-Lagrange equations (Ej 
and (Eq) correspond to classical mechanics in the potential 



Ur{x, y) = -]^x\a - -fy^f - ]^y\p - ^x^ . 



(28) 



When considering the optimal paths such as rc(i)j it is the potential Ur^ and not Vr, that is relevant. A plot of 
Ur{x, y) is shown in Fig. 5. 
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FIG. 5. Illustrative plot of the reduced potential UR{x,y) with a — 2, l3 — j = 1. 



We will, for concreteness, focus on paths which end in the positive x-valley; paths which end in the other valleys 
are treated in exactly the same way. The boundary conditions on these paths are a;(0) = y(0) = 0, x{T) — Xf and 
y{T) = yf, with yf <C Xf. We therefore simplify (p4) and (25) by keeping only terms which are linear in y, to obtain 
to lowest order, 



X = a X and 

j/=[/32-2(a + /3)7x2+7V] y. 



(29) 
(30) 



In our earlier treatment of this problem [y , we restricted attention to end-points lying on the x-axis: Xf = X and 
yf = (and correspondingly yf=Y and x/ = for paths ending in the positive y- valley) arguing that the relative 
contributions coming from paths ending on the axes should be approximately equal to the relative contributions 
coming from the paths ending anywhere in the valleys. In fact, as shown in Ref. jl|, this turns out to be an excellent 
approximation, but it has the deficiency that X is not determined. In this paper, we go beyond this approximation, 
which allows us to explore its validity, but also, as we will see, the precise value of X within this improved treatment 
does not enter; the result is independent of X, provided that X is large enough. The quantity j// has a different 
character to X: it will be integrated over at a later stage, when the probability flux through the valley is determined. 
The solution of ( p9| ) which satisfies the boundary conditions a;(0) = and x{T) = X is exactly the same as in Ref. 
0, namely. 



,{t) = X 



sinh at 
sinh aT 



(31) 
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The solution of ( |30| ) which satisfies the boundary conditions y(0) = and y{T) = yj (we reserve the notation Y for 
the final y value for the end-point of paths going along the y-valleys) is discussed in Appendix B. There it is shown 
that an excellent approximation to the solution of this equation (indistinguishable from the numerical solution) can 
be found by first solving the equation for small t and then for large t, and matching the two at some intermediate 
matching time i„i. Comparison with the numerical solution shows that the approximation remains good for a wide 
range of values of i„, from about O.IT to Q.6T or 0.7T. In order to be able to obtain a simple form for the solution, 



we have assumed that T is large, in the sense that e ""^ <C 1 and e ''"^ ^ 1. Using ( [B3| ) and (B7), we can write the 
explicit analytic from as 

r 2yf e-l^^ e''-^'/^" sinh/3i, < t < t^ 

yc{t) ^ i ~0(T^t) „ [-yX" 1^ „-2a(T-t)^] 4- ^ A- ^ rj. (32) 



y/ e ^^^ *^ exp 



^{l_e--(T-*)} 



tm <t<T 



For self-consistency we need to check that j/c(0 given by d32) is small, compared with Xc(f), if we are to justify the 
linearization procedure which led to ( p9|) and (pOJ). From (p2D it is straightforward to check that ycify has a maximum 
at i = t' , where t' = T, if X < X^^^ and T ~ t' ^ a^^ ln(X/Xmin), \i X > X,„in. As discussed in Sec. Ill, X,„in 
is the minimum value of x at which state selection can be said to have been completed, and therefore we ask that 
X > Xniin- The maximum value of ydt), 

increases rapidly with X and is already extremely large when X ~ bX^n^. It might therefore be tempting to argue 
that we need to take X to be much less than this value if the linearization procedure is to be valid. Since ( |33| ) has its 
least value when X — Xynin, this approximation would seem to be best for this choice of X, and in fact this was the 
choice made in our earlier work |l| . 

In fact, there are a variety of reasons why choosing X to be Xmin is not suitable. First, the argument above, namely 
that if X is too large then the linearization procedure is invalid, is in fact incorrect. We will find that the width of 
the distribution is so narrow in the y direction that the range of integration for yf need only be tiny — in fact only 
going out to values of y/ such that y? exp {7X^/0;} ~ 1. From ( |3^ ) we see that under these conditions yc{t) is indeed 
small. We will come back to this point again in Sees. VI and VII. 

The second reason is in fact more profound. Intuitively, we expect that for large enough values of X the probability 
flux through a given valley should be independent of X. This is due to the fact that in the regime of weak noise, 
once state selection has occurred there is no flux leakage out of a valley. This implies that we should not have to 
make any choice for X because our results should be independent of X as long as X is large enough. In fact, this 
is precisely what we find, namely that the probability flux through the x-valley tends to an asymptotic value as X 
tends to infinity. The remarkable cancellations that occur to produce an X independent fiux is a good indicator of 
the correctness of our calculational scheme. Again, we will come back to this point in Sees. VI and VII. 

Finally, throughout our calculation we use the approximation that T is large in the sense that e^"^ ^ 1 and 
g-pT ^ I This is apparently a problem, since we need the form of the distribution for all T in order to perform the 
time integral of the probability flux in the calculation of the total flux through the valley. The way out of this impasse 
is to take X sufficiently large that the probability current is essentially zero at small times, and only starts making 
appreciable contributions to the T integral when T is such that e~"'^ <C 1 and e~^^ <C 1. It is not clear how large X 
will have to be, but as discussed above, the probability flux through the x-valley tends to an asymptotic value as X 
grows, so we are assured of always finding a value of X for which the approximation that T is large is valid. 

At first sight, the large T approximation might appear to be problematic because it ignores the shorter state 
selection time scales a~^ and /3^^. However, this is not the case. The reason is that the integrals over time that 
are performed in the calculation of the action and Jacobian prefactor include these earlier times. In other words, the 
entire classical path, in particular the spatial region near the origin, is included in the calculation. 

After these rather technical asides, it is worthwhile summarizing what we have deduced about the e quat ion for the 
optimal path which leaves (0, 0) at t — and arrives at {X, yj) at time t = T. It is given by (|3^) and (^3) under the 
assumptions that (i) T is such that e^"-^ ^ 1 and e"''-^ ^ 1, (u) X is larger than a few ATmin, and (ui) tm and [T — tm) 
are large (we wifl later find that we require that e"^"*" < 1, e"^*^*'" < 1, e-2a(T-i,„) ^ i ^^^^ g-2/3(T-t,„) ^ -[.). 

We are now in a position to calculate the leading order contribution, PQ{X,yf,T), to the conditional probability 
distribution (^, by finding the action of the optimal path given by ( |3l| ) and (|3|). We begin by noting that, from 
(p6|), the classical paths are solutions oi x = —dU/dx and y = —dU/dy. Multiplying the first equation by i, the 
second one by y, adding them and integrating the result, gives ^(i^ + Vc) + U = E, a constant. Thus another form 
for the action of classical paths is 
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dt [xl + yl 



-ET 



\\y\l 



Details of the calculation are given in Appendix B. From ( |B20|) we find 



aX' 



1 



where 



bij 



f / a 



2a \jX 



file 



Sr. = ^ (cothaT -l) + ^yjSy., 



exp|2:^| r"" "dzz('5/")-ie-^(/3-az)2 + i [7^2-/3] 



(34) 



(35) 



(36) 



The two terms in (Kq) have different characters: as remarked earlier, we will eventually integrate over yf to ob- 
tain the probability flux through the valley, but the first term will remain, giving a leading order contribution of 
exp-(aXV4D)(cothaT- f) to (|2f]). 

When performing similar calculations to find the escape rate from one metastable state to another, the leading 
order result e^'^^^-^ gives a reasonable estimate, and the need to go to next order and to calculate the prefactor 
only arises if increased accuracy is required. The situation is very different in the present case of the decay from an 
unstable state. As mentioned already, we expect that the main contribution to the probability flux through the valley 
will occur when T ^ {2a)^^ ln{aX^/D). The conditional probability distribution ( pi] ) is peaked around a value of 
T of this order because of a balance between the leading (classical) term and the prefactor (a fluctuational) term. 
For this reason, the calculation of the prefactor in this problem, unlike in escape problems, is vital if the essential 
structure of the state selection probabilities is to be captured. We therefore turn to the calculation of this quantity. 



V. CALCULATION OF THE PREFACTOR 



In this section we will calculate the next to leading order contribution in (|2l|). From (|22| ) we see that it consists of 
two distinct contributions: from the Jacobian Jlf] evaluated at the optimal path, and from the fluctuations around 
the optimal path which give rise to the determinant of the second functional derivative of the action (|2^) . We begin 
by evaluation of this determinant. 

Starting from (|26|), and using the notation (|2^), we obtain 



Lire) 



' dt^ 



-ijxcyc [a + (3 - jxl 



-A-fXcyc [a + f3 - -fxl 



dt- 



+ {/3 - jxlY - 2ajxl 



(37) 



where, in line with our previous assumption, we have neglected terms which are 0(jj^) down on the terms shown in 
(p7|). We have omitted any time dependence in (|37| ) for the sake of clarity, but it should be remembered that the 
classical solutions Xc and yc are functions of t and the matrix is multiplied by an overall factor of d{t — t'). The 
diagonal entries are as in the yc — case |l) — only the off-diagonal entries are different. However, it is easy enough 
to see that the eigenvalues of (37) are even in y^ and since the matrix with yc — has no zero eigenvalues, the 
off-diagonal entries only provide 0{y'^) corrections to the eigenvalues of the matrix with yc = 0. Therefore neglecting 
these terms as before, we conclude that within the approximation we have adopted in this paper, we may set the 
off-diagonal entries in (|3^) to zero. From (g3|) and (^) we now find 

DetL(r;) = DetL^ Deti;, , (38) 

where L^ and Ly represent fiuctuations in x and y respectively: 
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d' 2 



(39) 
(40) 



To evaluate the determinants of the operators (p9) and (Eoh we make use of a well-known formula for such determinants 
p8[ . Let L be an operator of the type —cP/dt^^(t){t) and suppose that the eigenfunctions of L are required to vanish 
at the boundaries at t = and t = T, as in the case of interest to us here. Now let hi{t) and /i2(i) be two independent 
solutions of the homogeneous equation Lh ~ 0. Then 



Deti oc 



(41) 



/ti(0)/^2(T)-fe2(0)fei(T) 
/ii(0)/i2(0)-/i2(0)/ii(0) ■ 

The constant of proportionality in (|4^) may be omitted for the same reason the overall constant was omitted in (|2C 
the probability of ending up in the positive a;-valley will be normalized by the sum of the probabilities of ending up 
in the x- and j/-valleys. Furthermore, if /12 is the solution which vanishes at t = 0, then the formula for DetL involves 
only this solution: 



DetL = 



h2{T) 
/i2(0) 



The solution oi L^h = Q which satisfies h{Q) = is /i2(i) = Fsinhai, where i^ is a constant. Therefore 

sinh aT 
DetLi. = . 



(42) 



(43) 



To find DetLj, we need only note that the solutions of the homogeneous equation Lyh — are also the solutions of 
the classical equation (|30|). But the solution which vanishes at i = has already been found in Appendix B: it is 
given by (H) and (^), 



h2{t) 



Gsinhpt, 



f e''* exp 



-yX' ^-2a(T-t) 
2a 



< t < i,„ 

trr.<t<T 



where G is a constant. Since /i2(0) — /3G, use of ([4^ ) yields, 

1 



DetL,, 



e^'^ exp 



7^2 



2/3 ' """ [ 2q 
It only remains to calculate the Jacobian at the optimal path. From (lia) this is given by 

J{f,) = exp I i 1^ dt [(-a + jy^it)) + {-f3 + jxl{t))] \ . 

Neglecting the 0{y^) term as before, one finds that Jifc) — JxJy where 

1 



J. 


= exp j 


"2"^ 


Jy 


= exp<^ 


--\f3T 




w exp<^ 


-\PT 



Aa 
Aa 



sinh 2aT - 2aT 



2 sinh^ aT 



(44) 



(45) 



(46) 



(47) 



(48) 



since we are assuming that e^""^ <^ 1. Since both the determinants and the Jacobians factorize into contributions 
associated with the x coordinate and contributions associated with the y coordinate, we may summarize the results 
of this section so far as 



^x 



\/DetLx 



Ji, 



y/DctLy 



(2a)i/2g-aT^ 



(2/3)1/2 e-'^^ exp 



2^ 
2a 



(49) 
(50) 
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It is straightforward to check that these results agree with those of Ref. ||l| under the assumptions that e~°''^ <^ 1 
and e~^'^ <^ 1. While neglecting corrections to the leading behavior of the prefactor which are of the type e"^""^ or 
g-2/3T jg jjj jjj^g with this large T approximation, we will see later that the final integration over T will also provide a 
posteriori justification for the neglect of these terms. 



From (21), (p5|), (|4S|) and (pOj), we can now write down the conditional probability for the system to be at fj 



{X, yf) at time T, given that it was at the origin at time i = 0, as 
P(r>, r|0, 0) = (4a/3)^/' exp |-(« + /3)r+ ^ 



As we discussed in detail in Sees. II and III, when using the modified potential (S) the question is not what are 
the paths of least action from the origin to the minima of the potential, since those minima cease to exist when the 
limit (5, e ^ is taken. Instead, the question is what is the relative flux through one valley compared to the other. 
To calculate the flux, it is not the conditional probability distribution (|5l| ) that we need to know, but the probability 
current J^ — (^Jx^Jy)- This current is related to the conditional probability distribution through the Fokker-Planck 



equation |29 



dP 

— + divj^ = , (52) 

where 

J = -(Vy)F-i:»VP. (53) 

In order to find J from P given by (pl|), we first note from (^l]) that 

i:iVP= -fvp(°))p [1 + 0(1))] . (54) 



Therefore we only need to differentiate V and P*^°^ with respect to X and j/y in order to find J . Carrying this out 
we obtain 

Jx = (|X [coth(aT) - 1] + aX - -iXyf^ P (55) 

Jy^{SyVs^liVi--iVsX^)P. (56) 

In the next section we will calculate the flux through the positive x-valley. This will involve only the component of 
J normal to the line x = X. Therefore, the y-component of the current, jTy, will not contribute and need not be 
considered any further; the entire contribution to the flux will come from integrating Jx^ given by (p5|), over yj. By 
virtue of the y? term in the exponential in (|5l|) this will be a Gaussian integral. While we would naturally neglect 
the 0(yV) term in (|55|) in line with previous approximations in this section, we now see that it would in any case give 
a contribution of order D once the ?// integral is performed. This is a further justification for ignoring such terms. 
Finally, we have been neglecting e~^""^ type corrections in the prefactor, which means that we should replace cothaT 



by 1. These approximations lead to the result J^ — aXP and so we flnd from (51) 



Jx{X, yf,T) = 2a(a/3)i/2 X exp | -(a + /3)r + ^ | 



This is the key result of this section. We need now only use it to calculate the total flux through the positive oi-valley 
at X = X and compare it with the analogous quantity through the positive y-valley. This will give us the relative 
probability of the x-mode being selected. 
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VI. RESULTS 

The idea underlying the method we use to calculate the relative probability of one of the states being selected is 
most easily understood by first describing the analogous procedure used when carrying out Monte Carlo simulations. 
As was briefly alluded to in Sec. Ill, for each run the particle starts at the origin and subsequently follows the Langevin 
dynamics of the reduced problem until it reaches x = ±X in one of the x-valleys ot y = ±Y in one of the y-valleys. 
We used X — 5Xniin and Y — SFmin to be certain that state selection had occurred. At this point that particular 
run ceases. If the particle has ended up in a y-valley, for instance, 1 is added to the total number of runs which have 
selected the y-mode. Another run is initiated and the result of that is added to the totals. After a large number of 
runs the proportions selecting the x-mode and y-mode are used to calculate N^ and Ny, the relative probabilities of 
these states being selected. 

Let us focus on the runs ending up in the positive cc-valley, exactly as we have been doing in the analytic treatment 
in Sees. IV and V. For each of these runs the final value of j/ — called yj in the analytic treatment above — will 
vary. We would expect most of the runs to end near to the x-axis, with off-axis end-points becoming less and less 
common as we move away from the axis. This is reflected in the y/ dependence of jTx given by (|57|). Just as in the 
Monte Carlo simulation, where we add up all contributions with differing final y-coordinate at a: = X, so we need 
to integrate (p^) over all yj at x — X . Similarly, just as in the Monte Carlo simulation, where we add up all of the 
contributions, no matter how long they took to get to the end-point, so in the analytic treatment we have to integrate 
over all T to obtain the total flux. Therefore we need to calculate 



TAX)^ dT dyfJ.AX,yf,T), (58) 

Jo J -oo 

in order to calculate the state selection probabilities. 

While the above justification for ( p8| ) as the quantity which we need to calculate seems intuitively plausible, it is 
worthwhile formally proving this. We first need to specify what "ending up in the positive x- valley" means in terms 
of a mathematical expression. Since, in the limit T — > oo, trajectories will have entered one of the four valleys, we 
define the probability that it has entered the positive x-valley as 

dyf / dXP{X,yf,T\Q,0). (59) 

oo Jx 

Here the positive a;- valley is defined as the entire potential surface to the right of the line x — X. On the other hand 
from the continuity equation (p2) 



hm P{X,yf,T\6,0)-P{X,yf,0\0,0)= H ^P{X,yf,T\Q,Q)dT 
r^oo Jq ai 



divJdT. (60) 







Integrating (pG) over the "volume" {x > X, — cx) < y/ < oo}, using that fact that at T = the end-point is at the 
origin and making use of the divergence theorem, gives 

/oo /"OO /"OO /• 

dyf / dXP{X,yf,T\Q,0) = - dT J.dS . (61) 

-oo Jx Jo Js 

The "surface" integral on the right-hand side of (^) only gives a contribution at x — X, where dS is in the direction 
of the outward normal, i.e. in the negative x direction. Thus J .dS = —Jxdyf and so the required probability (|59|) is 
equal to ([5^). 

Substituting (p7|) into (pq) and performing the j/j integration yields 



:FAX) = 2a(a/3)i/2 X (^^) ''' ^^P {^} ^ ^T 



^<a+l3)T 



aX^ 
xexp<{ — — - (cothaT-1) J>[1 + 0(D)] . (62) 



As remarked in Sec. V, any 0{y^f) corrections in the prefactor give 0{D) corrections to (p2[), as do 0{yi) corrections 
to the action, which justifies their omission. The final integral over yf also provides justification for the linear 
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approximation. An examination of J'x{X,yf,T) shows that this function is effectively non-zero only for small values 
of |y/|, which get still smaller as X increases, so that unless \yf\ ^ exp {—■-fX^/2a}, the flux is essentially zero. Thus 
either the conditional probability distribution or the conditional probability current resemble a thin wafer centered 
on the X-axis when plotted at fixed T. Moreover, this wafer gets very much thinner with increasing X. This means 
that the limits on the y/ integral in (M) are effectively ± exp {—7X^/20;}, which is tiny for large X. It is more 
instructive to change variables from yf to y*f — yf exp {7X^/20;}, so that the integral now has the range (—1,1). 
But, as discussed in Sec. IV, this means that ( |33|) now reads 

which has a magnitude less than one for —1 < y*f < 1. Thus the linearization procedure is justified. In some sense, 
the variable y*f is more appropriate in this problem than yf itself. 

An examination of the integrand in (p2h shows that is has a maximum when T ~T* = (2a)^^ hi{a^X'^ /[{a + [3)D\), 
which is slightly different to the value mentioned in Sec. IV. This is because that value came from the estimates of Sec. 
Ill and Appendix A, which used the value for the maximum which is found in the analogous one-dimensional stochastic 
process in x{t). The more refined calculation we have been discussing here has an additional e'~"'^ contribution coming 
from fluctuations in y{t) and from Jy. It is the occurrence of this term in the integrand of (p2) which changes the 
value at which the maximum occurs by this very slight amount. From (p2) we can also check the conjecture made in 
Sec. IV, that if X is greater than a few Xn^in, then even the tails of the integrand will be at values of T which are large 
enough for the approximations made in this paper to be valid. A numerical investigation of the integrand shows that 
it has typically already fallen by two orders of magnitude from its maximum value for T satisfying e~""^, e~^'^ <C 1, 
when X is larger than a few X^i^. 

To perform the integration in ([6^ ) we (i) replace (cothaT — 1) by its large T form, 2e^'^°'^ (justified below), and 
(ii) change variables to ^ = {aX^ /2D)e~^^^ . Then the integral becomes 

1 / on \ ("+/3)/2a /-H 

where S = {aX'^)/2D. The integral is an incomplete gamma function: it equals r([a -I- l3\/2a) plus exponentially 
small corrections in D coming from the large, but finite, upper limit S p7[ . Neglecting these exponentially small 
terms we therefore obtain for the total flux through the positive a;-valley at X 

-'-) ^ <"«-'-^ (T)"^ -{.?} (1^)'"°"^" ^ m ^^^om . 

By using the approximation (cothaT — 1) r:! 2e~^""^, terms of order £)~ie~^'-"'^^-'""^, {n — 1, . . .) were neglected, but 
these are of the order of £)"^"+i and give contributions which are 0{D) down on the leading term (|65|). Similarly, 
any corrections to the prefactor of order e~"^ or e~^'^ would also give contributions which are down compared to 
(pq). Once again, we see that these flnal integrals justify approximations which we made earlier. 



The total flux given by (65) depends on X, which is still undetermined, apart from the requirement that it should 
be greater than a few X-min- However, if we investigate the X dependence of the expression in (£^), we find that as X 
increases from ^min, it first decreases until X equals two or three times X-aam where it reaches a constant value and 
subsequently remains at this value as X — j 00. Of course, this constancy of the flux is exactly what we would expect: 
once the system has selected the positive x-valley, it remains in that state and so the flux should be conserved within 
the valle y, i .e. be the same for all X. We will denote this constant value simply as Tx] it is given by the AT — > 00 
limit of (pa): 

/3/2q 
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The integral in (Bq) may be expressed in terms of gamma functions; it equals aj3T{P/a). Therefore 

P/2c 



Tx = 2D (2vr) V2 (hR) \ ^" / [1 + 0{D)] . (67) 
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The result ( pT] ) only depends on the four parameters Q;,/3, 7 and D, as we would expect. An exactly analogous 
calculation of the total flux through the y- valley at Y gives the same form but with a and f3 interchanged: 



J^y = 2D (27r) 
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a/20 r 
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(68) 



Throughout the analysis we have ignored constant factors which multiply both the results coming from paths ending 
up in the x- valley and those ending up in the y- valley. The reason given was that eventually we would normalize the 
probabilities of state selection and therefore common factors were irrelevant. We are now at the point where we can 
impose this normalization, but before doing so, let us put in the constant multiplying factor explicitly, and in doing 
so absorb the additional factors of 2D (2i^Y^'^ which appear in ( |67| ) and ( |6^ ) into it. We therefore the final forms for 
the total fluxes along the x- and j/-valleys as 
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where /C is the overall constant. 

The relative probability of flowing into an x-valley has already been introduced in Sec. III. In terms of the total 
fluxes it is given by 
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T + T ' 
with Ny = l-N,,. From Eqns. @), {^ and (0) we obtain: 
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where 
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Equations ^72) and (O) are the desired results. Note that the only place where the coupling constant 7 enters is 
multiplying the noise strength D. It is simple to see that this will always be the case since, if x and y are measured in 
units of X^in a-nd Ymin respectively in the reduced problem, the only place where 7 appears in the Langevin equation 
is multiplying D. Thus the effect of the interaction, specified by 7, is to renormalize the noise. This means that the 
probability of either the x-mode or of the y-mode being selected depends only on the three quantities a, (3 and 7Z?. 
It is clear that N^ and Ny given by (|7^), have the right behavior in various limits. First of all, if a = /3 they are 
equal. If a and (5 are not approximately the same, the relative magnitudes of N^ and Ny are largely governed by 
their dependence on D: if a > (3, then cr > 0, and N^ < 1 while Ny < D'^ ^ 1. A more detailed comparison with the 
results of Monte Carlo simulations is given in Fig. 6. 
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FIG. 6. Probability of flowing into an a;- valley as a function of a, with /3 : 
results and the full curve is our theoretical result for N^ given by (|72[). 



1,7 = 1 and D = 0.01. The dots show simulation 



The final result ([72|) is in good, but not perfect, agreement with simulations. We believe that the only restriction 
on our method is that D should be small, so that the method of steepest descent is appropriate. For D — 0.01, the 
neglect of 0{D) corrections to (p9) and (pOh should mean that our result is within approximately 1% of the simulation 
results. It is not clear from an examination of Fig. 6 whether the slight discrepancy is due to the 0{D) corrections 
or to additional contributions which have not been accounted for. Since the agreement is still not perfect for smaller 
D, we should take the latter explanation seriously and carry out a search to see if we can find such additional terms. 
A possible source is the following. We have conjectured that the solution of (pi) and (pSl) which has the least action 
can be obtained first by finding the Xc f£) solution with j/c — 0, and then linearizing about it. However, there are 
undoubtedly higher action solutions of ( |24| ) and (|2j) which have a non-trivial structure in x and y, that is, solutions 
which cannot be obtained as a simple expansion about the y — Q solution. In other applications of the steepest descent 
method, we would neglect these solutions as they would give contributions which were exponentially smaller than the 
contributions from the least action solutions, and thus be completely negligible for small D. However, in this problem, 
after integrating over T, the previously exponentially small contribution becomes a power law (see Eqn. (64)), or 
more precisely, D divided by the a quantity related to the action of the solution, all raised to a power. Although 
we do not know what power the higher action solutions will be raised to without further analysis, it is clear that 
solutions to (p4) and (Ea) with, for instance, twice the action of the solution considered in this paper, have at least 
the potential of changing our result by a not insignificant amount. Thus it appears that the neglect of higher action 
solutions, which is common in most problems involving optimal paths, may not be so well-founded in the situation 
which we are considering in this paper. 
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VII. CONCLUSIONS 

In this paper we have presented the motivation, justification and calculational details of a method for determining 
the probabiHties of various metastable states being selected for occupation, when an unstable state decays. Although, 
as we hope we have made clear, this is a widespread phenomenon, theoretical progress seems to have been hampered 
by a lack of sufficiently powerful tools with which to attack the problem. As a consequence, most studies were carried 
out some time ago on one-dimensional systems (or quasi-one-dimensional systems). Yet, these are actually the systems 
of least interest. The really interesting aspects of the problem arise in two or more dimensions, where the system 
can perform Brownian motion near an unstable state, and then select a final state through a mixture of randomness 
(noise) and dynamics (which will necessarily be non-linear in all non-trivial cases). It is this interplay of randomness 
and determinism, both essential to the process, which makes the problem so difficult. 

Our main purpose in this paper has been to show how the use of path-integral techniques, and especially the 
notion of optimal paths, can be used to analyze this technically hard problem. In order to stress that relatively naive 
approaches are not able to capture much of the subtlety of this problem, we carried through what seemed to us the 
most simple-minded approach to the problem in two dimensions, and showed that it was not able to give satisfactory 
results. Although the extent of the agreement with numerical work is always open to interpretation, there were two 
features which signified that the naive method was defective. Firstly, there was no systematic way to proceed with 
the calculation; one was faced with ad hoc choices which had to be made at several junctures in the calculation. 
Secondly, the crude approximation (^ to the already unsystematic result ( |l2[ ) , actually gives better agreement with 
simulations than (|lj) itself, which indicates the dubious nature of these results. It should also be borne in mind that 
the prediction that N^ = 1/2 when a = /3 is assured by symmetry and that in almost all schemes N^ ^ 1 as a 
gets very much larger than /3, simply because the noise will have very little effect in this case. If we add to this the 
expectation that Nx{a) will be a smooth function, we see that even the most simple minded approximation cannot 
be expected to be too far out. These observations indicate that we should focus more on the systematics (or lack 
thereof) of the calculation, rather than getting perfect agreement with simulation data. 

The nature of the optimal paths used in the treatment presented here is different from those in the more familiar 
problem of escape from one potential well to the other. In that situation it takes an exponentially large time 
(exp AF/D, where AV is the height of the barrier), to make the transition from one potential well to another, and 
this is reflected in the fact that the optimal paths are essentially of infinitely long duration (T -^ oo). In the problem 
studied in this paper we are focusing on phenomena on much shorter time scales: the decay of a state with no barrier 
to a metastable state. Any subsequent process (presumably a noise-activated escape to a different well, of the kind 
we have just discussed) is of no interest to us here. Actually, as explained in Sec. II, we do not even ask that the 
end point of the optimal path is a metastable state, but only that it lies somewhere in the valley which leads to the 
metastable state in the full problem. The optimal paths in these cases are of finite duration. To achieve this, the 
initial velocity of the particle in the mechanical analogy has to be non-zero at the origin. As T increases this initial 
velocity decreases, and eventually tends to zero if T tends to infinity. 

Many of the other details of the calculation involving optimal paths turn out to be remarkably subtle. For example, 
the use of the reduced potential simplifies the problem, but eliminating the actual minima means that we have to 
specify an arbitrary final value of x (if we are considering state-selection into the positive a;-valley) , which we denote 
by X. However, the probability of the positive x-valley being selected should not depend on X, as long as X has a 
value which lies in the valley. Reassuringly, this is what we find in our calculation, but it requires us to take into 
account the j/-dependence of the path for this "common-sense" condition to be observed. In our earlier treatment [y, 
we sought to give only the outline of our method, and made the assumption that the total probability flux through a 
particular valley was proportional to (a) its value on the axes, and (b) its maximum value, which occurs at T = T* 
— the time for which P{X, 0, 0|0, 0) is a maximum at a fixed X. Within this scheme there was no way of determining 
X and Y, and so we fixed them to our simulation results. The best values were X ^ A^min and Y -^ ymin, in line with 
our expectations. However, we have found here that it is possible to get a form for the y-component of the optimal 
path only if we assume that X is sufficiently large — which comes from the requirement that T must be sufficiently 
large — and so setting X « Xmin and Y « Fmin is not possible in the approach we have adopted here. Needless to 
say, we expect that if we were able to find the full solution for the optimum path, then we would not need to constrain 
A or y to be large: the expressions for the fluxes JF^ and Ty would be valid for all X > Xmin and for all Y > Fmin- 
Presumably the fact that !Fx decreases from its value at X = A„iin, and only becomes constant when X equals a few 
Amin, is not due to a violation of flux conservation (this would necessitate some flux leakage, and it is difficult to see 
where it would go to), but instead to the inapplicability of our calculational scheme for small X , which simply gives 
the incorrect expression for J-'^. It should be stressed that this is not a problem for our method: we simply take X 
and Y large enough so that J^r and Ty respectively, do tend to asymptotic values. These are the required values and 
are those which we use in our determination of N^ and Ny. 
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An even more subtle aspect of the calculation concerns the applicability of the linear approximation in determining 
the y-component of the optimal path. Although at first sight ydt) is not small for large X, it turns out that the range 
of final values of yc{t), denoted by y/, decreases so fast as X increases, that the small ydt) approximation is still 
valid. Put more intuitively: if T is large, the path has to take a wide detour to large y if it is to take sufficiently large 
time to reach {X,yf) (recall that T is fixed for these paths). But "large y" means large relative to yc{T) — y/, and, 
as discussed earlier, the current J^x{X, yf, T) is effectively non-zero only for tiny values of \yf\. So while it is true that 
on the scale of y/ the excursions of the optimal path spread deep into the plane, in fact the range of y/ is so minute, 
that the linearization assumption remains valid. The result that the probability is so concentrated about the axes is 
used in other places in the calculation. For example, in Sec. VI we argued that the "volume" integral in (pfl) could 



be replaced by the "surface" integral in (61), and then we proceeded to ignore some of the contributions from the 
"surfaces" in the later integral. We specified the volume in the case of the positive x-valley as the entire plane to the 
right of the line x = X. Actually, this is not quite correct, since we need to define a similar volume for the y-valleys, 
and these will overlap in the sectors where x and y are both large. We can be reasonably vague about how precisely 
to define these volumes and surfaces simply because Jx (and Jy) fall away so fast that they have an utterly negligible 
contribution however we define them. In any reasonable definition, the only contribution to the surface integral (|6l| ) 
will be on the line x = X , very close to the x-axis. 

We believe that one of the best indicators of the correctness of our approach is the cancellation of the various 
factors of X coming from several different sources: (a) the fluctuations about the yc (i)-solution, (b) the action of the 
yc(i)-solution (after integration over yj), (c) the action of the a:c(0"Solution (after integration over T), and (d) the 
definition of the current in terms of the probability distribution. What remains after these cancellations is essentially 
independent of X, as would intuitively be expected. As we have repeatedly stressed, many of these contributions 
come from assuming that T is large. This seems somewhat paradoxical, since we would expect that state selection 
is determined at earlier times. This can be understood if one realizes that the large X calculation determines the 
normalization of the flux in (p5h, or equivalently, the quantity (^(a, (3) in (|7^). To obtain a finite X ^ cx) limit, it is 
essential to include a non-trivial ydt) optimal path in the analysis. If the optimal path is taken to be ydt) — 0, then 
only the D'^ factor in ( [72| ) is determined. This was essentially the approach we adopted in our earlier paper [nl. While 
we believe that the treatment given here is a great improvement on that reported in Ref. |l| , there is no doubt room 
for further improvement. For instance, the linearization approximation seems very reasonable on physical grounds, 
but it would be useful to put it on a sounder mathematical footing. A deeper understanding of the origin of the small 
scale set by exp {—'^X'^/a) would also be valuable. 

The competition between new modes, when a previously stable mode becomes unstable, is responsible for much of 
the emergent order found in systems far from equilibrium. The continual branching to more complex structures that 
this entails pG], is not only governed by the dynamics of the process but also by the random fluctuations, or noise, 
generated by the large number of other degrees of freedom of the system not explicitly included in the description. 
The model which we have investigated in this paper is, as has already been emphasized, the simplest showing enough 
of the complicated features of this process of state selection to serve as an illustrative example of our method. In 
forthcoming papers we will show how the method can be applied to more complex examples, such as state selection in 
lasers and the study of population dynamics in a fluctuating environment. We believe that the ideas and techniques 
developed here will not only be applicable to these situations, but to many others where multiple states compete for 
occupation. 
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APPENDIX A: 

In this appendix, we explore the consequences of the elementary theory put forward at the end of Sec. III. We will 
obtain simple expressions for the probability that the particle ends up in an cc-valley, which are compared to Monte 
Carlo simulations in Sec. III. 

We begin with the first ingredient in this simple theory. Since 7 = we have only to solve the Langevin equations 
X — ax + Tjx{t) and y — f3y + Tjy{t). Such linear problems can always be exactly solved. Since x and y are linearly 
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related to rj^ and T]y, they are also Gaussian random variables. It is easy to show that {x{t)) = {y{t)) — and therefore 
the probability that the particle is at {x, y) at time t, given that it started at the origin at time f = is 



P(a;,y,t|0,0,0) =7Vexp 
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2{x^(t)) 2{yHt)) 
where A/" is a normalization constant. For later use we will also need 
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To quantify the second ingredient, first suppose that the coordinates of the particle {x, y) are in the first quadrant. 
If the angle it makes with the positive a;-axis, iaii^^ {y / x) is less than tan^^(ymin/A'ijiin), then according to our 
criterion, it goes into the a;-valley. Since the tan function is monotonic in the interval (0, 7r/2), and generalizing to 
other quadrants by using |a;| and |y| rather than x and y, the criterion becomes: 
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Thes e ingredients can now be put together. The question reduces to how often the particle is in the sector specified 
by (A3) at a given time t and how often it is in the sector specified by (A4) at this time. This probability is found by 
integrating ( |A1[ ) over all x and all y with the constraint that it satisfies either (A3) or (A4) depending on whether we 
want the probability of ending up in the x- valley or y- valley. Restricting ourselves to the positive quadrant, which we 
can obviously do on symmetry grounds, the probability of ending up in an x- valley, N^ and in a y- valley, Ny, are 
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where the factor 4 comes from including the contributions from the other quadrants and A = Y^i^/X^nin- In the N^ 
integral, change variables from y to z where y — Xxz at fixed x. The x integral can now be easily performed, giving 
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In a similar way 
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From (A2) we see that if a and (3 are not very different from each other, then K{t) will be neither very large nor very 
small, and we can ignore the integrals in ( A7) and ( A9) as a first approximation. Then determining J\f by asking that 
Nx + Ny ~ 1 and ignoring the factor 1 compared to the exponentials in (A2) yields 
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There are several simple argu ments wh ich give formulas resembling this. It is not difhcult to improve on this by 
including the integrals in (A7) and ([A9|). Once again determining A^ by normalization gives 
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2 2 

N^ = -tan"^[K(t)] ; Ny = - tan"^[K"\i)] , 

n IT 



(All) 



where we have used the resuh tan ^{w) + tan ^(1/w) = 7r/2. I f (a;^ (t)) » (y^(i)) or {y'^{t)) » (x^(i)), whic h 



happen even when a and /3 are not too different, we may replace ( All| ) by a simpler form which resembles ( A1C| ), but 
with constants multi plyin g the exponentials. 

The simple result ( AlO) does not depend on t he pa rameter 7, which represents the non-linear part of the potential, 
or on the noise strength, D. This is also true of ( All ) i f we ignor e the factor 1 compared to the exponentials in (|A2| ). 
To understand this, let us suppose that a> (3. From (AlC) or (All), we see that when t is small although there is 
a slight bias in favor of the a;- valley over the ?/- valley, it is not that marked. However, as t increases the asymmetry 
becomes more marked, until at large times the particle has only a very small chance of being in the y-valley. We need 
to simulate the role that the saddle points play in state selection by choosing the time to have the value tc, which is 
the time at which it is most likely that the x and y coordinates of the particle will have magnitudes X-^i-^ and Ymin 
respectively. This is implemented by asking that {x^{tc)) = X^^^ or {y'^{tc)) = ^min- Note that we have used the 
word "or", since only one relation is required to determine tc, and the two relations given give incompatible values for 
tc unless a — j3. Once again, as is common in many aspects of the naive approach to state selection which we have 
been summarizing in this appendix, there is a considerable degree of arbitrariness in the choice of ic- We will assume 
that for a > P, when the p artic le is most likely to go into the x- valley, tc is determined from the condition on (x^(ic) 



and vice- versa. Then from (A2) 
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which has the characteristic form — In I? as I? —» 0. 



If we substitute (A12) into (AlC) we obtain the probability of ending up in an x-valley. In terms of the quantities 
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we obtain 



N^ 



1 



1 + DP 



If we substitute (A12) into (All) we obtain the more general form 
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Equations ( A14) and ( |A15| ) are the principal results of this appendix. The two main points we wish to make is (i) 
their derivation is somewhat arbitrary — there are several other similar assumptions we could make which would give 
us slightly different formulas, and (ii) although, as is discussed in Sec. Ill, they show the right qualitative features, 
they are are not in good agreement with Monte Carlo simulations. Thus a more systematic and sophisticated theory 
is required. This is given in Sees. IV, V and VI of the paper. 



APPENDIX B: 



The purpose of this appendix is to explore the solution of ([3C|), where x{t) is given by ( pl| ) and where 2/(0) = and 
y{T) = yf. Throughout we will assume that T is large in the sense that e^""^ ^ 1 and e^^'^ <C 1, which will allow 
us to obtain a rather simple form for the solution. 



Small t solution. Let ( — X^cosech aT. Since e "^ <$^ 1, ( <$^ I. Therefore ( pOJ ) becomes 

y = y{(3^ - 2C7(a + 13) sinh^ at + 0(C^)} . 



We look for a solution of (Bl) as a power-series in C,: 

2/(t)=2/o(t)+yi(t)C + 0(C2). 
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Clearly the zeroth-order solution which satisfies y(0) = is yo{t) = Asinhpt, where A is an arbitrary constant. It is 
straightforward to determine yi{t), but since a short analysis shows that yi{t)( <C yo{t) for all t < T/2, we conclude 
that, it is sufficient to take the small t of the solution of (BOh — which we denote by ?/< — to be 
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y< (i) = A sinh /3t ; i < i„ 



(B3) 



Large t solution. In this regime both e ""^ ^ 1 and e "* ^ 1, and therefore x^(t) and a;f (i) may be approximated 
by X^e"^"^"^"*-* and X'^e^'*"*^^"*) respectively. Using these forms in (pOJ) and making the change of variables 
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we find that f{z) satisfies the equation 






This equation is easily solved to give 
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where A and B are arbitrary constants and where the upper linrit of the z-integral has been arbitrarily chosen to be 
the value that z takes on when t = T . 

We now have to match (86) to (B3) Sut t — tm- It will turn out that we do not have to specify t m pr ecisely, but 
let us assume for now that it has the value T /2. Then the corresponding value of z, z^ is seen from (B4) to be very 
small. Formally the matching procedure consists of equating both y{tm) and y(tm) for the solutions in the large and 
small time regimes. We will carry out this procedure below. However, let us first give a quick argument which gives 
us the correct final result. 

For t near tm, y<it) ^ Ae^^/2. Also from (B4), y>(i) = e^* e"^/^ fi^), and since z,n ^ 1, it follows that for t near 
tm, y>{t) ^ e'^* /(-z)- To get a smooth match, we require f{z) to be a constant, A/2, for values of z near z„i- Looking 
now at (p9), we see that the first term is constant, but the second term changes very quickly in the small z regime, 
since it diverges in the z ^ limit. We therefore take i3 = to get a smooth match, which in turn gives A = A/2. 
Therefore 
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Having matched (B3) and ( p37D near tm, we can now fix the constant A from the boundary condition y{T) = yf. 
Doing this gives us the form (|32|). 

Let us briefly indicate how we arrive at this result in a more systematic fashion. We can choose to either match 
{f{z),f'{z)} at z = Zm or {y{t),y{t)} at t = tm- They are related by 
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Now, by successive integration by parts, we can obtain an approximation to the integral in (B£) valid for small z. For 
our purposes we need retain only the first term: 
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Therefore to leading order 
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On the other hand y<{tm) == (A/2)e^*" and y<{tm) = {f3A/2)e^*"' . So from 
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Matching (|Bll|) and (|B12D gives 



A^^il + Oizra)] ; 6=^4"+^)/"[l + 0(z„ 



Substituting ( B13 ) into (B£) we obtain 
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Clearly for values of z that are not small, the second term is completely negligible and so />(z) ~ A/2 to a very good 
approximation. However, even for values of z near to Zm we can neglect this term — it gives a contribution which is 
0{zm)- Therefore, once again, we obtain (B7). 

In determining the form of y{t), we have implicitly assumed that the solution of least action is such that y{t) ^ 
for t 7^ 0. However, there are other solutions of the type y{t) — 0, for t < ti, (where ii < T) and a non-trivial solution 
of (pOJ), where x{t) is given by (31) and where y{ti) — and y(T) = yf, for t > ti. Intuitively, we would expect that 
the solution with ti — 0, which is the one found above, is the solution of least action, and all those with ti > have 
greater action. This is due to the fact that, especially for large ii, the path will have to sharply curve away from the 
X-axis if it is to satisfy the condition at t = T. We have checked this by numerically solving the differential equation 
( pO| ) for y{t) with different ii values, and found that the action monotonically increases with ti. 

Let us end by outlining the evaluation of the classical action (pj). To evaluate E we choose to take t = 0, since in 
this case U = 0. Also V\t^o = and V\t^T = -\aX'^ + \y)\l^~ P\- Therefore 







4 
1 



-aX' 



7ycL=o^ + >/[7^'-/3] 



(B15) 



The first bracket on the right-hand side of ( |B15D is easily evaluated from Xc{t) given by (|3l|). One finds 
(aX^/4)(cothQ;T — 1). The second bracket requires a little more calculation, since yc(i) has a different functional 
form depending on whether t is less than, or greater than, i^. Substituting the form for t < t„i first of all, gives for 
the second bracket: 



1 



A^(i 



A2/?2 



1 



dt yl -K -^ sinh 2/3t„ - -^ (T - t,„) + -2/^7^' - /?] 



(B16) 



Using (|32D and changing variables from t to z (given by (B4)) one finds 



^i>'^^<'>^T£(^) 



13/a 



^2(3T 



jX^/a 



dzz^^'°^^ 



{(3 - azf 



(B17) 



The range of the integral in ( B17 ) may be taken as (0, 7X^/0;) as long as we subtract out 

' dz z(^/")-i e"^ (/3 - azf = Oifiz^'' [1 -h 0{zm)\ ■ 



(B18) 



The u pper limit in the integral in (B17) is 0(1), and therefore the integral itself is also 0(1). Therefore th e con tribution 
( B18 ) is negligible comp ared to it, and so we may effectively replace the lower limit of the integral in ( B17 ) by zero. 
The third term in ( B16| ) being linear in (T — t,„) is negligible compared with the second term which is exponent ial in 
Ifitm- This term is in turn negligible compared to (B17), since g^^^^^^^*'") ^ 1. So the leading contribution to ( B17 ) 
is 



-/A^^'-^\^Wa\,X^ 



42 / ^ \/3/" j-iX'^ja 



(B19) 



Substituting the value of A determined from the boundary condition y{T) — yf, into (B19), the classical action (B15) 
becomes 
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1^2 1 

Sc^^- (cothaT - 1) + -y}bX^ - p] 



ill " exp^i^W dzz('5/")-ie-M/3-a^)'- (B20) 
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